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Direct numerical simulations of a single-step irreversible chemical reaction with 
non-premixed reactants in forced isotropic turbulence at R\ = 63, Da = 4.0, and 
Sc = 0.7 were made using 128 3 Fourier modes to obtain joint pdfs and other statis- 
tical information to parameterize and test a Fokker-Planck turbulent mixing model. 
Preliminary results indicate that the modeled gradient stretching term for an inert 
scalar is independent of the initial conditions of the scalar field. The conditional 
pdf of scalar gradient magnitudes is found to be a function of the scalar until the re- 
action is largely completed. Alignment of concentration gradients with local strain 
rate and other features of the flow were also investigated. 


1. Introduction 

Modern treatments of the theory of chemically reacting turbulent flows are often 
based on the probability density function (pdf) method, since in the pdf equations 
for the concentrations of the chemical species, the chemical reaction terms are closed 
in the statistical sense (O’Brien 1980, Pope, 1985). However, the mixing terms 
involving molecular diffusion are not closed, so statistical models are needed for 
these terms. The shortcomings of the commonly used coalescence-dispersion models 
and LMSE closures have been well-documented (Kosaly & Givi, 1987, Leonard & 
Hill, 1991), and more recent closures such as the mapping closure (Chen et al. 
1989, Pope 1991, Gao 1991) and the linear-eddy model (Kerstein 1991) are being 
investigated. 

In the present study, the Fokker-Planck (FP) closure is applied to the joint scalar- 
scalar gradient pdf for a two-species, single-step, irreversible chemical reaction 

A + B — ► Products (1) 

of non-premixed reactants in forced, homogeneous isotropic turbulence. The mass 
conservation equation for the concentration of reactant A ( 4>a ) is 


d<t>A d<t>A 

+ Ui 


dt 


dx i 


= D 


&*a 

dxidxi 


— Sa, 
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where Sa “ —k<f>A<f>B in the work described here (a similar equation described the 
concentration of species B ). In the current development, it is assumed that the 
scalar diffusivity D is the same for both reactants and the product of reaction and 
that all physical properties are constant, including the finite reaction rate constant, 
k. The joint pdf’s of the reactant concentrations and their gradients are used in the 
model discussed here to avoid problems common with closures based on pdfs of the 
reactant concentrations alone. 

The joint pdf equation for the scalars <f>A and <j> B and their gradients xpA and ips 
may be written 


dP(<j>A, 4 >b) 
dt 


9 - 3T- [SflP] 


d<f>B 


d<t>A 

D W a - D ^r [(^b\4>a,4>b)P] 
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— Molecular mixing terms 
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where P( ) is the probability density function of its arguments and the arguments 
to P on the right hand side of each equation are the same as those appearing on the 
left, Sa and S B are the reaction source terms (both equal to — k<f>A<pB in this work), 
dij = dui/dxj is the velocity derivative tensor, and the summation convention 
applies to repeated indices. Clearly, the reaction terms in the above equations 
are closed. In the traditional scalar pdf formulation involving only concentrations 
<p, the three mixing terms in (3) must be modeled. In the joint pdf formulation 
studied here, the modeling is postponed to the gradient pdf equation (4), wherein 
the three molecular mixing terms (not shown) must be modeled as well as the 
scalar gradient magnification terms, the last two terms in (4). In the development 
to follow, the FP model studied here is further simplified by considering only a 
passive progress variable <f> rather than both reactant concentrations and by treating 
the diffusion/reaction zones between the two reactants as locally one-dimensional. 
Among other things, this allows us to consider the magnitude of the scalar gradient 
\xp\ rather than the full gradient vector. 

In this work, we compare the results of stochastic simulations with results from 
direct numerical simulations (DNS) and sample the DNS results to evaluate various 
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quantities that appear in the pdf equations. Also, the computed fields were probed 
for physical insights suggested from previous simulations at lower Reynolds number 
(Leonard et al. 1988 and Leonard & Hill 1988, 1990, 1992). 

2. Approach 

2.1. Direct numerical simulations 

To provide data to check the FP model, a direct numerical simulation of station- 
ary, isotropic turbulence with chemically reacting scalars was carried out using the 
Rogallo (1981) method with 128 3 Fourier modes and low- wavenumber negative vis- 
cosity to provide the forcing. The turbulence was allowed to evolve until it reached 
statistical equilibrium, at which time scalar fields for the reactant concentrations 
were initialized and the simulations were continued. Two sets of reacting scalar 
initial conditions were used in the simulation. Case I was begun from “blob” initial 
conditions in which the two reactants are segregated into three-dimensional “blobs” 
with thin diffusion zones between them. The distribution of blobs was determined 
from a passive scalar field using a method similar to that used by Eswaran & Pope 
(1988). However, in this case, we follow Leonard & Hill (1991) and use a passive 
scalar field that has evolved with the turbulence so that the initial blobs are corre- 
lated with the velocity field. Case II uses ‘slab’ conditions, in which the reactants 
A and B are segregated into “slabs” with two planar ( x-z planes) diffusion zones 
between them in the periodic domain. In both cases, the overall (average) reac- 
tant concentrations were in stoichiometric proportions. The Damkohler number 
or dimensionless reaction rate coefficient was set at Da = kAoq 2 / e = 4.0 where 
q 2 = (u,-uj) and e is the dissipation rate of the turbulent kinetic energy 2 i/(e,^e, ; ) 

( eij is the strain rate tensor). The Schmidt number Sc was 0.7 for all species. 
The simulations were carried out until te/q 2 = 0.968. Figure 1 shows the reaction 
zones in the plane z — 0 at time te/q 2 = 0.968 for the two cases, showing the nearly 
isotropic scalar field for the blob conditions and remnants of the initial dual reaction 
zones in the slab case. 

Comprehensive diagnostics of the simulated fields were generated, including 
marginal, joint, and conditional pdf’s of the concentrations of the reactants and 
the conserved scalar 4> = <Pa ~ ^Bi the magnitudes of their gradients, velocity field 
properties such as the vorticity and dissipation, and various correlations. 

2.2. Gradient alignment analysis 

An analysis of the alignment of the reactant concentration gradients was carried 
out to provide theoretical support for stochastic models and closures that assume 
one- dimensionality or alignment of scalar gradients in non-premixed systems (this 
includes flamelet models, conditional moment closures, and the linear eddy model 
as well as the model developed here). The approach taken was to use (2) to ob- 
tain an expression for the evolution of the cosine of the angle n ab between scalar 
gradient vectors ip A and ipB , where hab = V’^iV’fli/dV’/illV’fll)- For cases in which 
Dfi/Dt = 0, a linear stability analysis was performed to determine the stability of 
this state. General results were obtained for arbitrary reaction rate functions S(</>), 
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FIGURE 1. Scalar 4>a4>b (reaction rate over k) in an x-y plane at te/q 2 = 0.96S 
for (a) blob initial conditions and ( b ) slab initial conditions. Contour increments 
are (a) 0.05 and ( b ) 0.1. Shaded areas indicate large values of the gradient 

amplification rate (ipAt^ij^Bj > 3.0 {ipAiCij^Bj))- 
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Figure 2. PDF of the cosine of the angle between and V’B at te/q 2 = 0.497 
for the blob case (fi ab = V , /»«V'Bi/(|V’^l|V’fl|))- 

including reversible reactions and the temperature-dependent kinetics. In addition, 
the alignment of reactant gradients with temperature gradient, the reaction product 
gradient, and the conserved scalar ( <f>A — 4>b) gradient was considered. The results 
pertinent to this study are summarized here: 

1. If IMb(< = 0) = -1 (gradients initially aligned and opposed in the non-premixed 
system) then hab{L) remains equal to —1 (aligned), independent of the reaction 
rate and of the presence of products of reaction (in the reversible case) including 
temperature, independent of the diffusivities Da and Db (which may differ), and 
independent of the strain rate e tJ , except as noted below. 

2. A stability analysis of the case described in (1) above shows that in nonisothermal 
systems, the reactant concentration gradients can become misaligned, depending 
on the Zeldovich number and on the direction of the temperature gradient with 
respect to ip a- 

3. If /xab( 0) — 1 (gradients initially misaligned) then the irreversible reaction (1) 

tends to align the gradients of <j>A and <f>B- 

4. As the reaction rate constant k in (2) becomes large, the reactants become seg- 
regated such that n — — 1 on the reaction surface and undefined elsewhere. 

5. The alignment of a reacting and non-reacting scalar, say <f> A and <f>, is preserved 
as in (1) above and is not influenced by the reaction rate even if gradients are 
initially misaligned. 

Thus, in the simple non-premixed reaction case considered here, the initial scalar 
gradients are aligned (ji a B = — 1), and remain aligned for all time. This theoretical 
result was confirmed in the direct numerical simulations by examining the pdf of 
HAb, which is approximately a delta functions at -1 (see figure 2). 

2.3 Fokker- Planck closure 

A Fokker-Planck (FP) molecular mixing closure was developed by Fox (1992a) to 
describe the evolution of the joint scalar, scalar gradient pdf in a system of reacting 
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one-dimensional, random-sized lamellae. Numerical (Fox, 1992b) and theoretical 
(Sokolov and Blumen, 1991a, 1991b) studies of diffusion in such systems have shown 
that the joint pdf evolves to a bivariate independent Gaussian pdf. However, if 
the scalar and scalar gradient are initially correlated, the correlation diminishes at 
a rate on the order of the scalar dissipation rate suggesting that the scalar and 
scalar gradient cannot be treated as independent random variables. Fox (1992b) 
has shown that the FP closure captures the form of the joint pdf and the decay 
rate of the correlation function for diffusion in random-sized lamellae in the absence 
of turbulent stretching, and suggested a modification to the closure to include the 
effect of the latter. In the following subsections, the application and extension of 
this model to the reacting system under consideration is presented. 

2.8.1 A single inert scalar 

In the following derivation for an inert scalar, the diffusion is assumed to be locally 
one-dimensional, so that only the magnitude of the scalar gradient is relevant. In 
§2.3.2, the alignment results of §2.2 will allow this treatment to be extended to the 
reacting multiple scalar case. For a scalar <f> and its gradient xp, the modified FP 
closure can be expressed in terms of a pair of stochastic differential equations: 

d<f> = k ip)dt + KB^(<f>y i p)dW <i> {t), (5) 

dtp = k 2 \p)dt 4- Cu*u*\pdt + KB^(<f>, %l))dW^{t)y (6) 

where A B and B $ are functions determined as in Fox (1992a), C w *LO*xp 
is the gradient stretching term suggested by Fox (1992b), lo* is the turbulence 
relaxation rate defined by Pope &: Chen (1990) (see (9) below), and 

k 2 = D(iP 2 )/(<t> 2 ) = 6D/ AJ. (7) 

The turbulence relaxation rate, a;*, is a random variable defined in terms of the 
(random) pseudo-dissipation rate, 


e*(x, t) — v 


du{ dui 
dxj dxj ’ 


(S) 


and the (nonrandom) turbulent kinetic energy, q 2 )/ 2 = (u,u,)/2, as 

U)*(x, t) = 2e*(x, t)/q 2 (t). (9) 

Pope and Chen (1990) have proposed a stochastic differential equation for u?* whose 
coefficients are independent of <f> and 0, and which yields a limiting log-normal pdf 
for 10 * . The gradient stretching constant, C w *, is assumed to be independent of the 
initial conditions. 

The FP model can be used to derive equations for the moments of the scalar 
and its gradient. In particular for an inert scalar in isotropic turbulence, the exact 
equations for the variance of the scalar is 


d(<{> 2 ) 


-2D(^ 2 ), 


dt 


( 10 ) 
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and, from the model, the scalar gradient variance is 

= -2CyD^- +2C w .(u;> 2 ), (11) 

at (<t> i ) 

where C y is a parameter in the definition of the functions Ay and By in (6). In 
the absence of turbulence (u>* = 0), the above equations are closed and constitute a 
two-equation model for the scalar energy and its dissipation rate. For this case, Fox 
(1992b) has found that Cy = 3 gives a good fit to the random-sized lamellae data 
and is required by the limiting form of (<p 2 )/(ip 2 ) predicted by theory (Sokolov and 
Blumen, 1991a). Note that if u>* and ip axe uncorrelated, or if u >* is nonrandom as is 
often assumed in pdf modeling studies, then (u>*ip 2 ) = (u*)(ip 2 ) and the equations 
are again closed. In particular, if (u>*) is time-independent, the long-time asymptotic 
behavior of the variances (characterized by constant (ip 2 )/{<p 2 )) can be determined 
as 

< 12 > 

Note that the scalar rms decreases exponentially in the limit of large t and the rate 
is independent of D. Other molecular mixing closures for the scalar pdf, such as 
the LMSE model (Pope, 1985) usually take 

^p = (is) 

with C $ = 2.0. 

i The FP closure discussed above extends standard scalar mixing models in two 

ways: (1) it models the scalar dissipation rate instead of assuming that is con- 
stant, and (2) it treats the turbulence relaxation rate as a random variable so that, 
for example, regions in the flow with large lj* will be correlated with regions of 
large scalar gradient and hence with increased mixing and reaction. Direct numeri- 
cal simulations indicate that these qualitative features are characteristic of turbulent 
reacting flows (see Leonard and Hill, 1991 and §3). 

2.3.2 Multiple reacting scalars 

! In the FP closure, multiple reactive scalars are handled by first considering an 

s inert scalar <f> with gradient ip, which are governed by the stochastic differential 

equations (5) and (6). Let <f> Q and tl? aj = 1,».,^» denote reactive scalars and 
their gradients, respectively, all with the same molecular diffusivity as and with 
] linearly dependent initial values; that is <j!> a (;r, y, z, t = 0) = a a <f>(x, y, z, 0) + b a , and 

if> Q (x y y, z, 0) = a a 4>{x,y,z, 0). As discussed by Pope (1985), the molecular mixing 
model to be developed below must be linear in the scalars <f> Q , so in the absence of 
source terms 

; ^ = F(^,^ 0 + Gi(f - (F(^, - (G,(^, ^ oi ), (14) 

: at 


i 
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where F find G, axe functions to be determined. Using the same arguments, it can be 
shown that a similar linearity property must hold for Moreover, when the source 
term is zero, <f> a (x,y, z,t) = a Q <j>(x,y,z,t) + b a and tp a (x,y,z,t) = a Q xJ?(x,y,z,t) 
for all t. The scalars thus remain correlated for all time, implying that the same 
two Wiener processes ( W $ and W^) that appear in (5) and (6) must also appear 
in the equations for 4>a and tp a (Fox, 1992a). Equivalently, the multiple reactive 
scalar model can be formulated in terms of d<j> and d\j> as follows. 

Assume that a local one-to-one differentiable mapping exists between <£ and </> Q , 
namely <f> Q = Differentiation then leads to an expression for the time- 

derivative of $ a : 


d<t> Q 

dt 






(15) 


Given <f> and Vn (15) is closed; however, it cannot be conveniently solved using Monte 
Carlo methods. 

It is interesting that the conditional moment closure (CMC) (Bilger 1991) employs 
a similar mapping equation: 


d* a 

dt 


= D(xl> 2 | 4>) 


d 2 * a 

d<\> 2 


+ S a ($i,...,$;v), 


(16) 


where D(\p 2 \<f>) is the conditional scalar dissipation rate of the inert scalar given <j>. 
The CMC mapping equation is closed given <j> and results in a joint pdf for <f) and <j> a 
with a 1-dimensional support (it falls on a curve in (<£, <^)-space). However, in the 
current formulation, the support will, in general, be 2-dimensional since each value 
of ip will result in a separate curve in (<f>, <^)-space. Mell et al. (1992) have solved 
the CMC mapping equation numerically for the reaction A + B — ► P with D(tfi 2 \4>) 
and the pdf of <f> taken directly from DNS, and found good agreement with the joint 
pdf of <j> and 4 M computed from the DNS data (the curve computed by CMC falls 
near the maximum of the joint pdf found by DNS). In addition, they found that the 
CMC results are insensitive to the functional form used for D(0 2 |<^), which suggests 
that the source term may dominate the diffusion term in the mapping equation. We 
therefore hope that the crude model used for this term below will be adequate. 

In order to apply a Monte Carlo method, d4> a is written in terms of its partial 
derivatives: 

d ^ = ^of d 4> + ^dt 1 (17) 

or substituting (15) 


d<t>c = -ZT # + w (-zir)dt + S a (4>u-,<f>N)dt. 


d<{> 


d<t> 2 


( 18 ) 


In (18) the diffusion term (premultiplied by D) is zero if the source term is zero or 
linear. Also, if the source term is such that 4> a is time invariant (e.g., an infinitely 
fast bimolecular reaction) the sum of the diffusion and source terms is zero. Oth- 
erwise it must be closed in terms of the <£, ip, 4 > a and i/> Q , and the closure must be 
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linear in <j> Q and xj> a as discussed above. The simplest closure hypothesis is to as- 
sume that this term is independent of the modeled variables, leading to the simplest 
model of the form (14): 


dK = *rf d4> ~ + s ^-^^ dt - 


(19) 


The terms involving di a /dtp are closed in terms of ip a (see below). 

To obtain a similar expression for dip ai , note that our assumed form for <p a implies 

that , d<t> Q di a , 

1pai = — = ST'I’i l 20 ) 


dx{ dtp 


The total derivative of ip a i is thus 




( 21 ) 


The time derivative term can be evaluated by differentiating (15) with respect to 
x i to obtain 


* dtd<f> 



ip 2 ipi + D 


d 2 i a d^ 2 

dtp 2 dii 



dS a dip 


( 22 ) 


which when substituted into (21) yields 


j, , ji , j , 

dip a =— —'Pi dtp + -st dtpi 


d<p 2 
+ D 


dtp 3 


dtp 

ip 2 ipi dt + D- 


Pi* ,,, „ . 


dtp 2 dxi 


j dSa dip , 


(23) 


Terms involving more than one derivative of 4>„ are not closed with respect to tp Q 
and i p Q . The first term is modeled as zero, which is exact if 5, is zero or linear and 
the diffusion terms are modeled as in (19) yielding 


dtpei = 


die > 
dtp 


dxpi - ( 


dia 

dtp 


d *Pi) + XI 

/» 


dS Q dip 

dtPp dtp 


ipi dt 


(24) 


The functional form of ip a , in (20) implies that all the scalar gradients are aligned 
with the gradient of the inert scalar <p , in agreement with the analysis in §2.2. 
Therefore, it is not necessary to treat ip a and ip as vectors. Without loss of generality 
we can let xp a = ipaiipi/^l (the projection of xp a onto the ip direction) and ip = \ip\. 
This is also consistent with the one-dimensional nature of the FP model for the 
inert scalar (equations (5) and (6)). It is also clear that di a /dtp = ipa/'P- With 
these simplifications the final model equations for the evolution of tp Q and ip a are 


dtp a 

dt 


Ipa d<P 
ip dt 


- + Sc,(<Pi,—,<Pn) 


(25) 
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te/q 2 


t*/q 2 


FIGURE 3. Evolution of the model constants (a) C^, and (6) C w . in the direct 
numerical simulation for blob initial conditions and slab initial condi- 

tions. 


d\p Q 

~df 


( — W dtp . dS Q . 

^ t b dt^° ^ %p dt ^ + ^ 


d<f>p 


(26) 


which in the absence of source terms are linear as required. 

Since the marginal pdfs of tp and \j> a are both symmetric about zero, the expected 
value of \j) obtained from (26) should be zero for all time. The FP closure for tf> Q 
is similar to the LMSE closure since, in the absence of source terms, d In \rp a \/dt = 
51n|0| /dt. But since d In \xj>\/dt is not constant, the two closures are not identical; 
in particular, with the FP closure the joint pdf of <f> and \j> evolves to a bivariate 
Gaussian (Fox, 1992a). 

The FP closure can be extended to scalars with nonequal molecular diffusivities 
by including a separate inert scalar with the same molecular diffusivity as each 
corresponding reactive scalar. The stochastic differential equations for the new 
inert scalars have the same form as Eqs. (5) and (6) except with k 2 modified to 
include the correct molecular diffusivity. The same Wiener processes (W$ and W$) 
must appear in each pair of stochastic differential equations as discussed by Fox 
(1992a). 


2-4 Evaluation of model constants 

The FP closure has two “universal” constants and C w * whose values can 
be determined using DNS data. The exact equations for an inert scalar and the 
magnitude of its gradient in a three-dimensional flow are 


D<f> 

Dt 


= DV 2 <p, 


(27) 


DxP 2 

~w 


2Dxpi 


<9 2 V>. 

dxj dxj 


2 faeijtltj, 


(28) 


where e,j is the strain rate tensor. By comparing the expected value of (28) with 
(11) derived from the FP model, the model constants can be evaluated as 
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FIGURE 4. Evolution of 1$$} evaluated from DNS for blob and 

slab initial conditions. 


and 


cv = - 


( 0 2 ) 2 


(30) 


The model constants, evaluated from the direct numerical simulations using (29) 
and (30), are shown as functions of time in figure 3. It is clear from figure 3a that 
the gradient mixing constant, C depends strongly on initial conditions. For the 
anisotropic slab initial conditions, Cj, decreases steadily with time after an initial 
jump. On the other hand, for the isotropic blob initial conditions, C\p ~ 6.7 for 
all time. This difference is most likely due the difference in integral length scales 
of the scalar fields in the two cases. Integral length scale effects have not yet been 
incorporated into the FP closure, but have been shown to have a significant effect 
on the scalar dissipation rate (Eswaran and Pope, 1988; Eswaran and O’Brien, 
1989; Kosaly, 1989; Jiang and O’Brien, 1991). For the slab case, scalar integral 
length scales are initially infinite in two directions and decrease as the turbulent 
advection creates a more isotropic field. For the blob case, the integral length 
scale is approximately constant for all time. In contrast, the gradient stretching 
constant, C w . , shown in figure 3b appears to be independent of initial conditions 
and approximately equal to 4.7 for all time. 

An expression relating C ^ and C^* can be found from the limiting value of 
A 2 = 6 (<j> 2 )/(rp 2 ). From Eqs. (10) and (11) for the FP model, the following relation 

is found for A 2 : 

^ = 12D(C* - 1) - 2 C u . ^^A 2 . (31) 


For statistically stationary w*, (31) has a single stable limiting solution: 


if) P(cy-i)(V’ 2 ) 

+' (i/> 2 ) C w .(u*rP*) ' 


( 32 ) 
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Thus, if the model has asymptotic behavior consistent with the DNS, we should 
find that 


DjCj, - l){fy 

Cw(u>*VW) 


(33) 


for large t. The left side of (33) is plotted as a function of time in figure 4 for both 
the slab and blob cases. For both cases this quantity does indeed approach 1. 

Mantel and Borghi (1991) have proposed a two-equation model for scalar en- 
ergy, scalar dissipation similar in form to (10) and (11) but with (w*)(^ 2 ) in place 
of In their model, the coefficients are defined in terms of a turbulence 

Reynolds number Re t = yj(q 2 /2)l t /v so that, for large Re t , C ^ = /? 0 \/Ke7/2, and 
Cv* = ocQs/Ret/2 , with a$ = 0.9 and /3q — 1.25 found by DNS (Borghi et a/., 1992). 
For Ret = 194 of the present DNS simulations, these expressions yield C ^ = 8.7 
and C u ,• =6.3. These estimates are both 30% larger than the values given above, 
implying that the ratio is similar. Since C ^ = 3 in the limit where Re t = 0 

(Fox, 1992b), the Reynolds number dependence embodied in these large- Re t rela- 
tions may not be valid for these relatively low Reynolds number DNS computations. 
(Although the Re t values for the DNS runs used to determine ao and (3q have not 
yet appeared in the literature (Borghi et al y 1992), they must be small due to the 
limitations of DNS.) This fact may explain some of the discrepancies between the 
two sets of values for the model constants. 


2.5 Application to the single-step reaction 

The FP closure described above has been used to generate joint pdfs of the 
reactant concentrations and their gradients for the 2-component, single-step reaction 
scheme (1) using a Monte Carlo simulation. For this reaction (25) and (26) yield: 


d<t>A = - (— d(j>) - k<j> A <l> B dt , (34) 

d<t>B = - k<t> A <f> B dt , (35) 

d$A = ~ - kMsdt - k<f> B il> A dt y (36) 

and 

dip B = - (^-difr) - k<f> A ipBdt - k<f> B xl) A dt . (37) 

The Monte Carlo algorithm uses fractional time-stepping to split the mixing and 
reaction steps into separate processes (Pope, 1985; Fox, 1992a). The expected 
values appearing in (34)-(37) are computed during the mixing step so that the 
mean values of the scalars and scalar gradients are constant during mixing. A 
constant turbulence relaxation rate, u ;* = (u?*) was used in (6). The resulting joint 
pdf’s are compared to the DNS results below. 
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FIGURE 5. Conditional PDF of ln(|vH) given <f> (P(ln(|0|) | <£)) for the direct nu- 
merical simulation with (a, 6) slab initial conditions and ( c,d ) blob initial conditions 
at (a,c) te/q 2 = 0.497 and ( b,d ) te/q 2 = 0.968. At these times <f> T ms Mms(* = 0) 
is (a) 0.846,(6) 0.645, (c) 0.426, and (d) 0.172. 


3. Results 

3.1 Statistics of the inert scalar field 

The DNS data have been used to compute a wide range of statistics involving the 
inert and reactive scalars and the magnitudes of their gradients. The marginal pdf 
P(<f>) for the inert scalar is nearly identical in shape, at a given rms value, to those 
reported in previous DNS studies. The marginal pdf of the log of the magnitude 
of the inert scalar gradient, P(ln |0|), approaches a nearly Gaussian form, but with 
a slightly longer negative tail, in agreement with the DNS results of Eswaran and 
Pope (1988). 
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FIGURE 6. PDF of V 2 <f> for the DNS with blob initial conditions at te/q 2 = 0.968. 


Of greater interest is the joint pdf of the inert or conserved scalar, 0, and its 
gradient, ip, which has been the subject of both theory (Bilger, 1976; Fox, 1992a, 
1992b; Gao and O’Brien, 1991; Meyers and O’Brien, 1981; Valiho and Dopazo, 1991) 
and experiments (Anselmet and Antonia, 1985). In most reactive mixing closures 
(e.g., the flamelet model), <f> and ip are assumed to be independent so that the joint 
pdf is separable, P(<f>,tp\t) = t). In order to check for independence 

using the DNS data, the conditional pdf of ifi given <^, defined by P(%p\<p\t) = 
P(<f>, 0; t)/P(<f>] t)y has been computed. Note that if <j> and xp are independent, 
then P{xp\<p\t) — P(ip\t) and will thus be independent of <f> . Examples of the 
computed conditional pdf are shown in figure 5. From figure 5a (<^rms/<£rms(^ = 
0) = 0.846) it can be seen that near the start of the molecular mixing process, 
the scalar and scalar gradient are correlated. The correlation decays slowly so 
that in figure 5b (^ r ms/^rms(0) = 0.645) the conditional pdf continues to show a 
important ^-dependence. Not until the molecular mixing process is farther along 
(^rms/^rms(0) = 0.172) as shown in figure 5d does the conditional pdf appear to be 
independent of <f>. 

Another interesting pdf is that of the Laplacian of <f> shown in figure 6 on a log- 
linear plot for <£rms/^rms(0) = 0.172. As is clear from the form of the pdf, it is 
non-Gaussian with nearly exponential tails over 4 decades. 

3.2 Statistics for inert scalar mixing 

Statistics involving the scalar or scalar gradient and various turbulence quanti- 
ties have been computed using the DNS data. For example, the scalar gradient- 
turbulence relaxation rate correlation function, defined by 


2 ) 

(»' )( 4 > 2 ) 


(38) 


was found to be approximately time independent with values of 0.06 for the blob 
case and 0.15 for the slab case, indicating that u;* and have a slight tendency 
to be simultaneously larger than their mean value. This tendency can be seen 
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ln(u>*) 

Figure 7. Conditional pdf of ln(|^|) given ln(u> *) for the DNS with slab initial 
conditions at te/q 2 = 0.968. 



ln(u;*) 

FIGURE 8. Conditional pdf of <f> given ln(u;*) for the DNS with blob initial 
conditions at te/q 2 = 0.968. 

more clearly by examining the conditional pdf of ln(|^|) given ln(u;*) shown in 
figure 7 (slab case with <Arms/<£rms(0) = 0.645). From this figure it can be seen 
that the conditional pdf has a nearly constant shape but shifts upward as ln(u;*) 
increases. This behavior is consistent with the FP closure (6) wherein u* appears 
as a stretching (positive) term in the drift coefficient. 

Similar conclusions can be drawn from the conditional pdf of <f> given ln(u;*) shown 
in figure 8 for the blob case at <Arms/^rms(0) = 0.172. There it can be seen that 
larger values of ln(u>*) lead to smaller conditional variances for <j>. This is consistent 
with the model equations in that large ln(u;*) leads to large gradients and hence 
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FIGURE 9. Conditional expectation of <f>A<f>B (reaction rate over k ) given en- 

strophy (cj 2 ) and dissipation (e 2 ) for (a) slab initial conditions and ( b ) blob 

initial conditions at te/q 2 — 0.968. 

faster scalar dissipation. 

Finally, as noted by Pope and Chen (1990), the DNS simulations confirm that 
the log of the pseudo-dissipation rate of the turbulence, lne* is more nearly Gaus- 
sian than is the log of the true turbulence dissipation rate, In e. For example, the 
skewness and flatness of lne* are —0.06 and 3.05, respectively, compared to —0.29 
and 3.24 for In e at te/q 2 = 0.968 

3.3 Statistics of reacting scalars 

Some statistical quantities evaluated in previous simulations in decaying turbu- 
lence at lower Reynolds numbers (Leonard et ai 1988, Leonard & Hill 1988, 1990, 
1992) and in a similar study for a non-reacting system (Nomura Sz Elgobashi 1992) 
were examined in order to determine the extent to which the present system exhibits 
the same physical behavior. For example, pdf’s of the cosine between the directions 
of the reactant scalar gradients and the eigenvectors of the strain rate tensor, and 
plots of the eigenvectors superposed on reaction rate contours, show that there is 
considerable tendency for the most compressive eigenvector to align with the scalar 
gradients and to lie across the reaction zone. Furthermore, there is a similar but 
less pronounced tendency for the intermediate strain rate eigenvector to lie tangent 
to the reaction zones and isoscalar surfaces. 

Figures 9 and 10 show the effect of certain kinematic quantities on the reaction 
rate, and vice versa , at f=0.92. In figure 9, for reaction rate conditioned on levels of 
strain and enstrophy, (<f>A<f>B I e 2 } and (<j>A<f>B | u; 2 /2) where e 2 =e^e,j, it is seen that 
strain has a marked effect on reaction rate, but the effect of vorticity is considerably 
less. The converse plot, figure 10 for strain and enstrophy conditioned on reaction 
rate, confirms the previous observation of Leonard & Hill (1990) that conditional 
averages of e 2 and u> 2 /2 are near their volume averages and each other, except for 
the regions of most intense reaction rate where the straining is very high and the 
enstrophy is appreciably less than the volume averaged value. 

Regions where the gradient amplification term, (ipAi^ij^Bj)^ is greater than 
3.0 (tpAi^ij^Bj) are shaded in figure lb. Clearly the largest values of this term 
are associated with peak values in the reaction rate, supporting earlier claims by 
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FIGURE 10. Conditional expectation of enstrophy (u; 2 ) and dissipation 

(e 2 ) given (reaction rate over k) for (a) slab initial conditions and (6) blob 

initial conditions at tejq 2 = 0.968. 


Leonard &; Hill (1990, 1992) and underlining the importance of the modeled gradient 
amplification term in the FP model. Although not shown, the gradient amplifica- 
tion term seldom takes on negative values and since xj>A and x \>q tend to be aligned 
and opposing in the reaction zone, the compressive part of e,j must dominate this 
term as expected. 

Finally, various scalar gradient-strain rate correlation coefficients important in 
mixing studies were evaluated. One such quantity, {x^Ai^ij^Bj) / {^Ai^Bi) %/e 1 2 , ap- 
proaches the value —0.45 for the slab case and —0.40 for the blob case; the same 
values are obtained for the conserved or inert scalar <j> in these two cases. These 
values differ somewhat from the values —0.56 and —0.45 (—0.52 and —0.43 for the 
nonreacting scalars) found in decaying turbulence by Leonard & Hill (1990) and 
the value -0.5 found by Kerr (1985) for a nonreacting scalar in forced turbulence. 

The joint pdf of the reacting scalars, 4>a and <f>B have also been computed using 
the DNS data and can be compared to the joint pdf found from the FP closure. For 
example, the joint pdf at <£rms/Vrms(0) = 0.426 is shown in figure 11 and that found 
using the FP closure for the same value of <£rms/<Arms(0) in figure 12a. It can be seen 
that, despite the closure approximations needed to derive (19) and (24), the general 
shape of the pdf predicted by the FP closure corresponds closely to that found by 
DNS. In particular, the width of the curved region of significant probability is about 
the same in the DNS and the model. 1 

The comparisons between the pdfs of the modeled and DNS gradients (figures 
ll(b-d) and 12(b-d)) are not nearly as good, though the modeled gradients do have 
the correct order of magnitude. The strange bimodal structure of the gradient- 
gradient pdf (figure 12d) is presumably caused by one of the modeling assumptions 
used to derive equation (24). 


1 This can also be compared to the CMC model which predicts no scatter about the curve (Riley, 

1992 ). 






FIGURE 11. Joint pdf’s of reactant concentrations and gradients from the DNS 
for blob initial conditions at te/q 2 = 0.497. 


4* Conclusions 

Direct numerical simulations of a single-step chemical reaction between non- 
premixed reactants in forced isotropic turbulence were made for both “slab” and 
“blob” initial scalar reactant configurations. As found in previous simulations at 
lower Reynolds number, the amplification of concentration gradients in the reaction 
zone by the strain field was seen to be an important feature of these flows, in that 
regions of large local reaction rate are coincident with regions of large values of the 
gradient amplification factor. 

An analysis of the alignment of various scalar gradients with each other pro- 
vides some justification for treating the mixing process as locally one-dimensional 
as assumed in the Fokker- Planck model studied here and other closures. 

Comparisons were made between predictions of the FP closure and results of 




_4H 1 1 1 1 -4 T ' 

0.0 0.5 1.0 1.5 2.0 -4 -2 0 2 

<j>A H\^A\) 

FIGURE 12. Joint pdf’s of reactant concentrations and gradients from the FP 
model for the same conditions as in figure 11. 

turbulence simulations. The closure’s treatment of gradient stretching as a bilinear 
term in the model equation is generally supported by the DNS data. For example, 
the gradient stretching constant was found to be independent of initial conditions, 
and the DNS results for the joint pdf of the scalar gradient and the turbulence 
relaxation rate were found to be consistent with the model. Likewise, the closure’s 
prediction for the joint pdf of the reactive scalars is very similar in shape to the DNS 
result. However, it was also found that for the non-isotropic initial scalar field the 
gradient mixing constant appearing in the closure is not constant as assumed, and 
that the closure’s prediction for the form of the joint reactive scalar gradient pdf 
differs significantly from the DNS result. The former can most likely be accounted 
for in the closure by incorporating scalar integral length scale information, and the 
latter by modifying the closure assumptions used in deriving (24) from (23). In any 
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case, since it can be easily incorporated into existing Monte-Carlo simulation codes 
(Pope, 1985), the formulation of the FP closure in terms of a stochastic process 
offers a significant computational advantage over other closures that require the 
solution of a reaction- diffusion equation. 
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